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Abstract 

We introduce a new class of dynamic point process models with simple and 
intuitive dynamics that are based on the Voronoi tessellations generated by the 
processes. Under broad conditions, these processes prove to be ergodic and produce, 
on stabilisation, a wide range of clustering patterns. In the paper, we present results 
of simulation studies of three statistical measures (Thiel's redundancy, van Lieshout 
and Baddeley's J-function and the empirical distribution of the Voronoi nearest 
neighbours' numbers) for inference on these models from the clustering behaviour in 
the stationary regime. In particular, we make comparisons with the area-interaction 
processes of Baddeley and van Lieshout. 



1 Introduction 

In many diverse fields such as biology, geography, business and telecommunications, spa- 
tial point configurations evolve according to criteria dependent on a zone of influence in 
some general sense. A natural way of formalising this concept is based on using Voronoi 
tessellations (see e.g. P). 

Let (M, d) be a metric space and > the initial number of points in our process. 
The associated configuration space ^ consists of all A^-point subsets of M: ^ := {x C 
M : card(a;) = A^} ('multiple points' will be a.s. impossible in our models). For an 
X E ^ , the Voronoi cell of Xj relative to x is defined as 

■■= \y eM : d{y,Xi) = mind(|/, Xj)|. 

The set = {C^^ : Xj G a;} is called the Voronoi tessellation generated by x. 

We define a Voronoi point process as a discrete-time Markov processes {xn}n>o with 
values in which evolves as follows: at each step, a point is chosen from the current 
configuration x^ according to a probability rule determined by and removed from the 
configuration, and at the same time a new point is added at a random location, according 
to a fixed probability measure /i on M. Our initial interest interest in such dynamics 
was prompted by its relevance to some of the real-life processes in the above-mentioned 
application areas. In addition, simulations showed that processes constructed according 
to this model display very interesting forms of clustering behaviour, and constructing 
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models of point processes producing desired clustering patterns is of substantial interest 
both theoretically and practically, (see e.g. 

Despite the fact that it is notoriously difficult to obtain analytic results for processes 
of such a nature, we were able to prove ergodicity in a number of interesting cases. To 
formulate the results, we first need to further specify the 'culling rule' in the model. 

We assume that a non- negative 'selection function' S is given on the set of all possible 
Voronoi cells in M. Then, given that x = {xi, . . . ,xn} is the current configuration of the 
point process, in the next step a random point xj E x is chosen to be removed, with the 
distribution of the index J given by 

The function S can be based on different properties of a Voronoi cell. In this note, we 
consider only two of them: the 'volume' of the cell and the number of its edges. More 
precisely, we introduce the following two classes of Voronoi point processes: 

(A) The volume-based Voronoi point process, or i;-process, on M = §^ or M = [0, 1]^. 
Let A be the 'volume measure' on M (length on area on [0, 1]^). We assume that 
the value of the selection function S{C^.) is determined by the volume of the cell C^.: 
for a function : M"*" i— > M"*", one puts S{C^.) := Sy{X{C^.)). If Sy is increasing, then 
points with Voronoi cells of large volume are more likely to be culled, and so the selection 
favours points with small cells, i.e. points restricted by 'close neighbours'. A decreasing 
Sy favours points with large cells. Functions of the form Sy{u) = a G M, produce 
scale-independent models (in this case the dynamics of the process will clearly be invariant 
under scale transformations of M). 

(B) The neighbour-based Voronoi point process, or n-process, on M = [0, 1]^. The set 
x[xj] of Voronoi-nearest neighbours of a point xj G a; is defined as the collection of those 
generators G a; whose Voronoi cells share an edge with C!^. : 

x[xj] := |a;fe G a; : k ^ j, card(C,"^, H C,"^) > l}. 

We assume that, for a function Sn '■ {!,..., A^} ^ M^, we have S{C^.) = Sn{caTd{x[xj])). 
In both cases (A) and (B), we assume the placement probability fi = X. 

Theorem 1. (i) The neighbour-based Voronoi point process with a positive selection func- 
tion Sn '■ {!,..., A^} ^ is Harris ergodic. 

(ii) The volume-based Voronoi point process with a selection function Sy : (0, A(M)] 
M"*", such that both Sy and 1/Sy are bounded on closed subsets, is Harris ergodic. 

Recall that Harris ergodicity entails convergence to the stationary distribution in total 
variation (see e.g. p. 560 in jH] or p. 154 in For the proof of the theorem, see [S]. 

This result justifies the approach taken in the present note, which is devoted to simu- 
lation studies of point patterns emerging in stationary regime in the evolution of Voronoi 
point processes. More specifically, we are interested in the performance of three relatively 
simple statistics one could employ to characterise the resulting stationary distributions. 
These are Thiel's redundancy measure (which is essentially the relative entropy for the 
set of cell volumes) and, in the two-dimensional case, the distribution of the number of 
Voronoi nearest neighbours for a random cell and the J-function of Baddeley and van 



Lieshout (which provides a comparison of the environment of a 'random point' of the con- 
figuration with that of a 'random point' of the underlying space over a range of scales). 
We demonstrate that the n- and f-processes produce (for different choices of parame- 
ters) quite different clustering patterns, and that the combination of the above three 
measures work reasonably well in distinguishing between different Voronoi processes, and 
also between Voronoi processes and the area-interaction processes of Baddeley and Van 
Lieshout [3]. 

2 Simulation Studies 

2.1 One-dimensional t'-processes 

The simple model of scale- free process on M = with Sv{u) = m° demonstrates that 
very interesting dynamics and point patterns arise even in one- dimension. This includes 
a 'phase change' observed when varying the value of a. Figure [T] depicts side-by-side 
realisations of three different w-processes of this type, each having the same number of 
points = 128, but different a values. The base of each rectangle represents the circle 
opened out into a line segment by a cut. 



Figure 1 : Evolution of t)-processes on circle. Each of the three processes (with the selection functions 
S^{u) = with (a) a = -1.0, (b) a = 0.5, (c) a = 1.5, resp.) was run with TV = 128 points for T = 4096 
steps (the vertical axis represents time). 

We see a dramatic phase shift in behaviour between (b) and (c) that can be located 
more precisely at a = 1; its sharpness was found to increase with A^. When a > 1, a 
cluster forms whose stability also increases with A^. When a < 1, we observe a degree of 
clustering in M that varies with a. 

When a = 0, the points are uniformly distributed, so the associated Voronoi cell 
volumes have pretty nearly a Gamma distribution with shape parameter 2 (as cell-width 
is half the sum of two consecutive uniform spacings). For a < 1, the histogram of 




cell volumes stabilises after a large number of steps to a result well-fitted by a Gamma 
distribution with a shape parameter varying inversely with a. This suggests that the 
maximum likelihood estimator for the Gamma shape parameter for the cell-volume data 
could be a suitable statistic for inference on a. Furthermore, under the assumption 
of Gamma distribution, the shape parameter can equivalently be estimated from the 
entropy — and the latter can also be used without the above assumption, namely in the 
form of Thiel's redundancy measure which is introduced as follows. 

For a fixed configuration x, let pj := X{C^.)/X{M), 1 < n < N, he the probability 
that a random uniformly distributed point in M lies in C^.. Then Thiel's redundancy 
measure R*{x) of x is defined [Z,i2j as the entropy of the distribution {pj}jLi relative to 
the uniform one on {!,..., N}: 

N 

R*{x) := InN + y^^pj\npj. 

Simulations show that this statistic works quite well over a reasonably large range of a 
values: there is a very small variance in its values when > 10^, and a strong enough 
dependence on a to use R* for reliable estimation of the parameter. Figure|21 displays 
(connected) average values of 25 independent realizations for each of the values a = 
-2 + 0.3A;, < < 10, of R*{xt) after a large number of steps T (cf. Fig.Ej). 
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Figure 2: Empirical Thiel's redundancy measure for 7;-processes on §^ with S.i,{u) — u"-, N — 10"^. The 
dotted lines show the 95% probability intervals. 

The behaviour of the statistic R*{xt) as T increases could also be used to estimate 
the rate of convergence to stationarity. Figure|21 (left) gives an indication of that rate for 
a range of a values. It shows that the rate is quite high for a < 0.5, whereas when a 
approaches one, it takes the process much longer to settle in a stable regime, and also that 
the oscillations of {R*{xt)}t>o are higher for those values of a. At the threshold value 
a = 1, the sequence oscillates in quite diverse ranges of values demonstrating metastable 
beahviour, see Fig. El (right). 

2.2 Two dimensions: edge effects 

The problem of 'edge effects' (where the region boundary proximity can affect the shape/ size 
of Voronoi cells) in simulations of Voronoi tessellations in (a rectangle) M C is com- 
monly dealt with by treating the underlying space as a torus or by considering only 
those parts of M that lie within a window significantly distant from the edge dM of M. 
This approaches is satisfactory for 'static' point processes [9 . In the case of the Voronoi 
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Figure 3: Evolution of {R*{xt)}t>o for ?;-processes on S^: N = 10^, Sy{u) = u". Left: a = -1; -0.53; 
-0.05; 0.19; 0.45; 0.66; 0.9, on the time interval T < lOiV. Right: a = 1, on the time interval T < lO^A^. 

processes, edge effects may propagate in the course of the evolution, which presents an 
interesting problem by itself. 

Let a; C M = [0, 1]^ be finite, A C M^. The Voronoi tessellation induces a nearest- 
neighbour (NN) distance d^j^{xj,A) between Xj G x and A, defined as the length of the 
shortest 'path' xj = . . . , x{i) G x with the properties that the Voronoi cells of x{m) 
and x{m + 1) have a common edge, 1 < m < i, and C^{i) H A 7^ 0. We can gauge the 
significance of edge effects by observing the changes in R* when we restrict attention to 
Xj G X with d%j^{xj, dM) > m for a fixed m > 0. Edge effects were assessed by a two-way 
ANOVA using NN distance as one factor, and the selection function as the other. In the 
case of the w-process with S^{u) = m", we found that edge effects were only significant in 
the layer of cells adjacent to the boundary, without any significant effect from a. For the 
n-processes, significant edge effects were detected only for the 'anti-few' and 'anti-many' 
selection functions (see subsection 12 . 31 below) . with no significant propagation beyond the 
NN-depth of two. 

To minimise edge effects in our results, the statistics were generally computed from 
the set of cells with d^j^{xj, dM) > 3. 

2.3 Two-dimensional n-processes 

In this family of Voronoi processes, the evolution is driven by one of the most natural 
local characteristics: the number of neighbours of the process points, i.e. the number of 
the edges of their Voronoi cells. For non-degenerate random configurations in M^, the 
vertices of the Voronoi tessellation a.s. terminate three edges, so by Euler's theorem the 
mean number of edges for a cell is six. There can be cells with any number of edges 
m > 3, although in our simulations cells with m > 13 were rare, with relative frequency 

We studied the n-processes on M = [0, 1]^ with the following selection functions: 
(i) Sn{n) = n ('vanilla', in which cells with a large number of neighbours are more likely 
to be culled); (ii) Snin) = ('anti-many', which more severely penalises cells with large 
numbers of neighbours); (iii) Sn{n) = (n — 2)~^ ('anti-few', which does just the opposite); 
(iv) Sn{n) = (0.1 + \n — 6|)^^ ('anti-6', which penalises cells with six or close to six 
neighbours); (v) Sn{n) = |n — 6p ('pro-6', which does the opposite); and 'sharp filters' 
focussing on cells with a given number neighbours: (vi) S'n(5) = 5000 and S{n) = 1 if 
n 7^ 5 ('anti-5'); and (vii) S'„(5) = 1 and Sn{n) = 5000, if n 7^ 5 ('pro-5'), and in addition 
the 'anti-' and 'pro-' selectors for four and seven. 



The statistics used in the study included: 

(a) the empirical probability mass function (EPMF) for the number of Voronoi NN's, 

(b) Thiel's redundancy measure, and also 

(c) the J-function of van Lieshout and Baddeley 

The J-function J(r) compares the 'environment' of a 'typical generator' of the process 
with that of a 'random point' in the underlying space, and is defined as the ratio of 
the probabilities that a disk of radius r centred at the given point is empty of (other) 
generators of the process. 

Formally, let s be a stationary isotropic point process in M^. Set 

B{r) := {x eR^ : \x\ <r}, F{r) := P{x n B{r) 0), G (r) := {x n B (r) 0) , 

where is the reduced Palm distribution of x (the conditional distribution of ic\{0} 
given that there is a point at the origin). Then the J-function is defined as 



J(r) 



1 - G{r) 
1 - F(r) 



for all r > such that F(r) < 1. 



(2) 



In the case of a Poisson point process of constant intensity, or complete spatial randomness 
(CSR), clearly J(r) = 1, see e.g. If there is a tendency towards clustering obesrvable 
at scale r, then J(r) < 1, while J(r) > 1 if there is a tendency toward regular spacing 
and r is of the order of the average distance between neighbours. 

To estimate J{r), one uses empirical distribution functions F{r) and G{r): 



J{r) 



G{r) 



1 - F(r) 



This is only reliably computable when F{r) < 0.85. Moreover, simulations showed that 
the values of J have relatively large variances when x 10^, so we averaged the J values 
over 25-75 independent draws of the process, and performed a weighted cubic regression 
on the resulting points (with weights equal to the reciprocals of the standard deviations; 
higher order regressions add little to the goodness of fit). 

The results of our study for a selection of the above-listed n-processes are summarised 
in Fig. m (statistics (c), (a)) and Table[T] (statistic (b)). 



n-process 




s.d.(i?*) 


s.e.(:R') 


vanilla 


0.1450 


0.0068 


0.0014 


anti-many 


0.1548 


0.0073 


0.0015 


anti-few 


0.1074 


0.0047 


0.0009 


pro-6 


0.1615 


0.0083 


0.0012 


anti-6 


0.1524 


0.0059 


0.0008 


pro-5 


0.1479 


0.0067 


0.0013 


anti-5 


0.1255 


0.0063 


0.0013 



Table 1: R* values for various rt-processes, with standard deviations and standard errors. The values 
Tt are averages of 25 independent draws of R* (N = 2000, after T = 12N steps). For CSR, R* = 0.135. 

The J-function curves provide the most nuanced 'picture' of the distribution patterns, 
but these curves were obtained by averaging over a number of observations of highly 
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Figure 4: Comparison of n-processes {N = 2000, after T = 12N steps). Left: smoothed plots of In J(r). 
Right: Proportions of cells with given numbers of neighbours (given at the right-hand end). Standard 
error is approx. 1% of the given value in each case. 

variable statistics and are not easily applied in practice. The statistic R* provides a crude 
but fairly sensitive measure of deviation from CSR. The NN EPMF is probably the most 
effective statistic for the n-processes, but it doesn't distinguish well between CSR and 
'anti-many' models (and also between 'anti-few' and 'pro-6'), whereas the J-curves for 
these models are very different. The J-curves for these and for the 'sharp filters' on four 
and seven neighboured cells reveal the interesting fact that selection functions favouring 
cells with m > 6 (m < 6) Voronoi neighbours produced more (less, resp.) uniform [than 
CSR] configurations. The 'anti-few' and 'anti-5' processes produce very similar J-curves, 
but are distinguished by their effects on R*: 'anti-few' gives R* = 0.107, while 'anti-5' has 
R* = 0.126, which are significantly different from each other as well as being significantly 
below the CSR value of 0.135. The low value for R* corresponds to less variablity in cell 
areas, so this is consistent with the upward J-curves. Notice also the unusual curvature of 
the 'anti-6' curve, a result which was confirmed by increasing the number of independent 
draws in this case to 75. In all these cases, the R* and J(r) statistics gave consistent 
results, i.e. upward curves corresponded to R* < 0.135. 

2.4 Two-dimensional t'-processes 

We study here the f -processes on M = [0, 1]^ with the scale-free selection functions 
Sv{u) = u°'. FigureEl shows a few examples of point patterns after a considerable period 
of evolution. There is a continuous gradation of patterns, from almost regularly spaced 
ones for a < through increasing levels of clustering until we reach the threshold level 
at a = 1. For a > 1, the patterns become increasingly concentrated in a small region, 
usually in a corner of the square, as the number of steps increases. 

For a < 1 we can employ the same statistics as for the n-process. As in the one- 
dimensional case, Thiel's redundancy measure R* serves well for inferring the value of a, 
and also provides a good indicator of the stabilisation of the process, see Fig.lHl From 
Figs.iniElit appears that the two-dimensional processes settle into equilibrium at roughly 
the same rate as the one-dimensional processes. The In J(r) curves for a G (—3, 1) form 
a fan as expected, see Fig.[7|(left). More detailed analysis reveals surprising (statistically 
significant) changes in curvature that occur for a between —0.3 and —0.2 and between 0.5 
and 0.6: it appears that a subtle phase change does indeed occur around these values. In 




Fig. [7| (right) we have plotted the NN EPMF as a function of a. Note that, as a increases 
(and so the amount of clustering), the proportions of cells with different numbers of edges 
tend to come together; in other words, the cell geometries become more variegated. 




Figure 6: Thiel's redundancy measure for w-processes on [0, 1]^ (N = 2000, Sv(u) = -u"). Left: Evolution 
of {R*{xt)\t<12N for a = -3; -2.2; -l;-0.2; 0.2; 0.6; 1. Right: R* vs a (after T/N = 12 steps). 



3 Comparison with the area- interact ion process 



The area-interaction point process (AIPP) of Baddeley and van Lieshout [3] is a model 
which successfully produces a range of different clusterings. Kendall presented a prac- 



Figure 7: Graphical summary of results for u-processes on [0, 1]^ with Sy(u) = u" for different a values 
{N — 2000, T = 8N). Left: Smoothed curves for In j(r) (numbers of independent draws range from 
9 (a < —1.0) to 59 {a = 0.5,0.6)). Note the rapid change in the curve as a approaches 1. Right: 
Proportions of cells with given numbers of neighbours (given at the right-hand end) as functions of a. 

tical method for 'perfect simulation' of this process. This is summarised in Ambler 
from whose detailed account we derived our computer program for simulations. 

Let M = [0, 1]^ and ^ be the space of all finite configurations of points in M. The 
AIPP is specified by the density of its distribution with respect to that of the unit rate 
Poisson process on M. Let A be Lebesgue measure on M, G := B{p) for a fixed p > 0, 
and X (B G := {z : z = x + g,x E x, g & G} for x G The above density is given by 

p{x) = C/^^^'^dC^j^-AC^eG)^ 

where /5,7 > are parameters, G a normalising constant. Note that the Poisson process 
on M with constant rate /? has density ^card(a;)_ rpj^^ parameter 7 defines the interactive 
component of the process: when < 7 < 1, configurations with high values of X{x © G) 
for a given card(a;) are favoured, and so the interaction is described as repulsive, while 
7 > 1 favours x with low X{x Q) G) and so constitutes the attractive case. 

The expected number of points in M is a complicated function of 7, /? and p. We 
chose p = 0.01, 7 = 7^° with 71 G [0.3, 1.5], and adjusted j3 so that the number of points 
produced in the simulations of the AIPP was close to 2000. This is because the AIPP 
doesn't re-scale in a simple way, so comparison with the Voronoi u-process, in particular 
in regard to the J-function, requires the samples to have roughly the same number of 
points. 

Figures ISHTUl show the values of our three statistics for the AIPP. We observe a con- 
tinuous range of degrees of clustering on either side of CSR (when In J(r) = 1). It is 
interesting to compare these with those derived from the Voronoi f-process. In particu- 
lar, we can match results which produce close values for R*. We have done this for the 
cases of the AIPP with 71 = 1.5 and Voronoi f -process with a = 0.5 (for both processes, 
R* ^ 0.2). The right-hand plots in Figs.lHl HUl show that there is a significant differ- 
ence in the J-function curves, but that the Voronoi NN EPMF's do not yield significant 
difference. 

In summary, the u- processes and the AIPP's produce a range of clusterings depending 
on a continuous parameter. The difference in the resulting point patterns can be detected 
by the combination of our three statistics (a)-(c). 
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Figure 8: Plot of R* vs 71 for the AIPP as described above, with cubic regression Hne. Standard 
deviations of R* are close to 0.005. 
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Figure 9: Left: Plots of smoothed In j(r) curves for the AIPP's with various values of 71, showing a 
range of attractive/repulsive effects. Right: The graph of In J(r) for the AIPP with 71 = 1.5 is compared 
to that for the scale-free Voronoi t;-process with a = 0.5 (these two processes produce almost identical 
R* values and NN EPMF's). 
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Figure 10: Left: Proportions of cells with given numbers of neighbours (given at the right-hand end) 
for the AIPP's, plotted vs 71. Right: Comparison of the NN EPMF's for the AIPP with 71 = 1.5, the 
Voronoi u-process with a = 0.5, and CSR. Standard errors are close to 0.005. 
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